Physiological responses to drought stress of three pine species and comparative transcriptome analysis of Pinus yunnanensis var. pygmaea

Drought stress can significantly affect plant growth, development, and yield. Fewer comparative studies have been conducted between different species of pines, particularly involving Pinus yunnanensis var. pygmaea (P. pygmaea). In this study, the physiological indices, photosynthetic pigment and related antioxidant enzyme changes in needles from P. pygmaea, P. elliottii and P. massoniana under drought at 0, 7, 14, 21, 28 and 35 d, as well as 7 days after rehydration, were measured. The PacBio single-molecule real-time (SMRT) and Illumina RNA sequencing were used to uncover the gene expression differences in P. pygmaea under drought and rehydration conditions. The results showed that the total antioxidant capacity (TAOC) of P. pygmaea was significantly higher than P. massoniana and P. elliottii. TAOC showed a continuous increase trend across all species. Soluble sugar (SS), starch content and non-structural carbohydrate (NSC) of all three pines displayed a "W" pattern, declining initially, increasing, and then decreasing again. P. pygmaea exhibits stronger drought tolerance and greater recovery ability under prolonged drought conditions. Through the PacBio SMRT-seq, a total of 50,979 high-quality transcripts were generated, and 6,521 SSR and 5,561 long non-coding RNAs (LncRNAs) were identified. A total of 2310, 1849, 5271, 5947, 7710, and 6854 differentially expressed genes (DEGs) were identified compared to the control (Pp0D) in six pair-wise comparisons of treatment versus control. bHLH, NAC, ERF, MYB_related, C3H transcription factors (TFs) play an important role in drought tolerance of P. pygmaea. KEGG enrichment analysis and Gene set enrichment analysis (GSEA) analysis showed that P. pygmaea may respond to drought by enhancing metabolic processes such as ABA signaling pathway, alpha-linolenic acid. Weighted gene co-expression network analysis (WGCNA) revealed GST, CAT, LEC14B, SEC23 were associated with antioxidant enzyme activity and TAOC. This study provides a basis for further research on drought tolerance differences among coniferous species. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-024-10205-5.


Introduction
Pinus massoniana Lamb.(P.massoniana), mainly distributed from 21°41′N to 33°56′N and 102°10′E to 123°14′E [1], plays a vital role in soil and water conservation, restoration, and pulpwood production.The broad distribution and diverse uses make P. massoniana pivotal for sustainable forestry production and ecological environment construction [2].Compared with Pinus elliottii (P.elliottii), P. massoniana grows slower in the early stage (10-15 years) and faster in the later stage (> 15 years) [3].P. elliottii is an important timber and resin-producing species grown worldwide [4].In China, due to its strong adaptability, high resin yield, low resin crystallization rate, and high turpentine content, P. elliottii has become one of the main tree species for resin-tapping [5].P. yunnanensis suffers from genetic degradation, with high proportion of twisted trunk, leading to the emergence of mutated variants such as the Pinus yunnanensis var.pygmaea (P.pygmaea) [6].P. pygmaea maintains its dwarfism even growing in a suitable habitat [7].Compared to angiosperms, gymnosperms are characterized by higher vulnerability to subsequent droughts and more evident lagged drought effects ("legacy effects") [8].The abovementioned three pine species have a widespread distribution and significant economic value in China.However, there is limited research on the response to adversities and stress among different pine tree species.
Frequent natural abiotic stresses such as low temperature [9], drought [10], and low phosphorus [11] make us realize that the strength of tolerance is a key factor that affects and limits the normal growth into forests, and directly determines the survival rate and stand retention rate after afforestation.Due to climate change, drought will become increasingly frequent, severe, and longer [12,13].Drought significantly reduces the photosynthetic rate, the accumulation of biomass, and inhibits plant growth [14,15].Plants respond to drought stress through morphological and structural changes, expression of drought-resistant genes, hormone synthesis, and osmotic regulation substances [16].Under continuous drought conditions, trees close stomata or reduce leaf area to reduce water lose [17].For instance, P. palustris responded to drought treatment with greater stomatal control of plant water loss [18].Plant hormones, especially ABA, are considered to be closely associated with drought responses.Stomatal closure is controlled by both active (abscisic acid(ABA)-mediated closure) and passive (hydraulic closure) regulation [19].ABA-dependent stomatal closure is among the most rapid plant reactions to water stress [20].ABA accumulates during drought stress and rapidly decreases upon rehydration [8].Longterm water shortage induces the formation of stress memory that maintained a higher level of ABA accumulation in trees grown in more arid long-term conditions [21].Osmoprotectants such as total soluble sugars and proteins were early accumulated in primed plants during the stress [22].Meanwhile, antioxidant protection systems perform an essential function during drought stress [23].Many cell components are damaged by reactive oxygen species (ROS) that are generated as a result of decreases in photosynthesis [24].Antioxidant enzymes, such as catalase (CAT), peroxidase (POD) and superoxide dismutase (SOD) play a central role in plant adaptation to environmental changes [25].Osmolarity-related compounds and antioxidant enzyme activity will increase under drought stress [26].The high level of antioxidative enzyme activity is important for plants to tolerate stress for keeping good growth state and forming good quality [27].Plants can adapt to prolonged drought by remodeling transcriptomes [28].The recovery ability after rehydration is important for the successful adaptation of plants to arid environments [29].Recovery of P. massoniana seedlings is more likely to be inhibited when plants experience increasing drought stress [30].
Third-generation Pacific BioSciences (PacBio) singlemolecule real-time (SMRT) sequencing does not need to interrupt RNA fragments, and direct reverse transcription can be used to obtain full-length cDNA [31,32].Full-length transcripts that do not require assembly can be used to detect genes with AS events, APA events, lncRNAs and fusion transcripts [33].Pacbio SMRT sequence is widely used to obtain long reads and assemble a high quality reference transcriptome.PacBio-SMRT sequencing can provide a high-quality reference transcriptome for non-genomic species.RNA-seq is one of the most important and practical methods for exploring conifers genes [34,35].RNA-seq helps to provide differential gene expression information under different stress conditions, time periods, and genotypes [36].Resilience strategies differed among tree provenances: wet forests showed higher growth resistance to drought, while dry forests presented faster growth recovery [37].The response to stress pressure varies inconsistently among different species.Pterocarya stenoptera and P. elliottii seedlings developed different adaptive strategies in response to continuous flooding [38].Among P. elliottii, P. palustris and P. taeda, P. elliottii is the least droughttolerant [39].There has been extensive research on the drought response of P. massoniana.P. pygmaea has narrower niche and more adaption for extreme ecological conditions [40].However, fewer comparative studies have been conducted between different species of pines, particularly involving P. pygmaea.The physiological and molecular mechanisms of P. pygmaea are unclear.
In this study, the needles of three pines (P.pygmaea, P. elliottii, P. massoniana) were used as the research material to measure the antioxidant enzyme activity and physiological responses after drought and rehydration.By combining PacBio SMRT seq and RNA-seq, the molecular responses of P. pygmaea were studied to identify key genes and pathways.This study provides a foundation for breeding drought-resistant pine trees, and comprehending the physiological and molecular mechanisms underlying drought response.

Physiological responses of three pine species under drought and rehydration
Physiological indicators were measured in the three pine species (P.pygmaea, P. elliottii, P. massoniana) under prolonged drought stress.As the drought time lengthened, total antioxidant capacity (TAOC) showed a continuous increase across all species.At the Re7D, P. pygmaea had significantly higher TAOC than P. massoniana and P. elliottii (Fig. 1A).The soluble protein (SP) content of P. elliottii exhibited an initial decrease followed by an increase, with levels significantly higher than those of P. massoniana and P. pygmaea at several points between 14 and 35 days of drought stress (Fig. 1B).The soluble sugar (SS) content (Fig. 1C), starch content (Fig. 1D), and Nonstructural carbohydrate (NSC) (Fig. 1E) of all three pines displayed a "W" pattern.Initially, they decreased under prolonged drought stress, then increased, and finally decreased again before recovering after rehydration.After rehydration, P. pygmaea's POD enzyme activity was higher than the other two species (Fig. 1F), and malondialdehyde (MDA) content was lower than P. elliottii and P. massoniana at Re7D (Fig. 1G).During drought, the content of Chlorophyll a (Chla) and Chlorophyll b (Chlb) decreased continuously (Fig, 1H&I).After the drought treatments on three pine species, compared to the control (0d), TAOC showed an increasing trend, while chla and chlb exhibited a decreasing trend (Fig. 1J).Principal component analysis (PCA) (Fig. 1K) indicated that Re7D samples of the three pines were relatively similar, while the control (CK) samples of the pines were closed to each other.; K PCA analysis.Note: In (A-I), 7-35D represented drought treatment after 7 D, 14D, 21D, 28D, and 35D.Re7D represented 7 D after rehydration; error bars indicated the standard deviation, the number of seedlings measured for each species was n = 3; the left Y-axis represented value of physiological traitors, "CK" group of each species was set as a reference to compare the mean by the Student's t-Test, "ns" means no difference, "*" P ≤ 0.05;"**" P ≤ 0.01;"***" P ≤ 0.001;"****" P ≤ 0.0001.A least significant difference test (LSD) was used at a probability level of 0.05 to verify the significance between species at sample point

Screening and functional annotation of differentially expressed genes (DEGs)
For Illumina RNA-seq, a total of 45 Gb raw reads were obtained, the Q20 and Q30 scores for all samples of clean data aboved 99%, and the average content of GC was 44.54%, indicating good sequencing quality.Then, the clean reads were successfully mapped back to the full-length transcripts of PacBio SMRT.PCA analysis showed that samples within the same group were close to each other, indicating good repeatability (Fig. 3A).Based on the DEGs screening criteria, DEGs among all pairwise combinations were identified (Fig. 3B).Among them, compared with the PpCK (Pp0D), the numbers of DEGs between Pp14D and Pp35D were 1,849 (1082 (58.52%) up-regulated and 767 (41.48%) down-regulated) and 7,710 (3,486 (45.21%) up-regulated genes and 4,224 (54.79%)).There were 541 DEGs commonly shared between all groups compared with the PpCK (Fig. 3C), and 499 DEGs commonly shared between all groups compared with the Pp35D (Fig. 3D).Cluster analysis of the DEGs shared by PpCK revealed that genes within C1 (2,708 gene) and C2 (1,564 gene) were highly expressed in PpCK, but lowly expressed in all other samples which were related to drought stress.Genes within C6 (1,960 gene) and C7 (1,218 gene) were highly expressed in PpRe7D and maybe related to rehydration (Fig. 3E).

Identification of co-expression module
All filtered DEGs were retained for WGCNA analysis.The analysis identified 15 co-expression modules (expect for grey module) (Fig. 5A), of which bisque4 module was Fig. 4 STEM analysis of the gene expression pattern and GSEA analysis.Note: In A-H, each subfigure was a significant profile.The upper part of the subfigure showed the number and corresponding trend of the modules.Black line represented the profile trendline.Each line represented the relative expression trend of a gene.In I, GSEA of KEGG pathway between Pp35D and Pp0D; In J, GSEA of the "abscisic acid-activated signaling pathway" among the four groups; In K, GSEA of KEGG pathway between PpRe7D and Pp35D; In L, GSEA of the "plant-type cell wall" among the four groups positive correlation with POD (correlation coefficient (r) = 0.93, p = 1.3 × 10 -9 ), while salmon and white module was negative correlation with POD, with r of -0.55 (p = 0.0093), -0.51 (p = 0.019), respectively.The brown module was negative correlation with TAOC (r = -0.75,p = 1 × 10 -4 ).The plum2 module was positive correlation with TAOC (r = 0.7, p = 0.00037).The genes within the bisque4 module showed relatively lower expression during the term of rehydration (Fig. 5B).The genes within the plum2 module showed a higher expression at Pp28D and PpRe7D (Fig. 5C).The genes within the brown module showed a lower expression after Pp14D (Fig. 5D), while the darkgreen module genes showed a highest expression at Pp28D (Fig. 5E).

Validation of RNA-seq expression analysis
To validate the accuracy of RNA-seq data, seven genes were randomly selected for the quantitative real-time PCR (qRT-PCR) analysis.The qRT-PCR results showed that the expression patterns of genes were similar to their results in the RNA-seq (Fig. 7A-G).A strong positive correlation (r = 0.91, R 2 = 0.73) was obtained by a linear regression analysis (Fig. 7H), suggesting that the transcriptome data was reliable.

Discussion
In the present study, we investigated the response of P. pygmaea, P. elliottii, and P. massoniana seedlings to drought and subsequent recovery.We measured changes in the antioxidant system and oxidative stress levels, and evaluated the physiological plasticity of these species following rehydration.
Fig. 5 Weighted gene co-expression network analysis of the genes A Correlated heatmap of the adjacency of modules; B-E The eigengene expression of bisque4, plum2, brown, darkgreen module, respectively.Note: In a, each row represented a module, the color and number of each cell represented the correlation coefficient between modules and traits, the top number in the cell represented the correlation coefficient, the bottom number represented the p-value, the module name was shown on the y-axis, the color scale on the right shows module-trait correlation from −1 (blue) to 1 (red) Fig. 6 Partial of the co-expression network of module bwtween structure genes and transcription factors.A bisque4 module; B plum2 module.Note: In A, the top three genes which filled with different colors represented structural genes, with 'Pp' abbreviating 'D1_transcript_' .The TF families were shown below.In B, 'Pp' abbreviated 'D1_transcript_' , and the genes which filled with different colors in the middle represented structural genes.The arc encompasses genes of the same class, with TF connected on the periphery.In both A and B, the points within each gene family were randomly distributed inside the circle, with the radius of the circle proportional to the number of distribution points.The outermost points were calculated and connected using the chull function of grDevices R Package.The connections between structural genes were mapped in 'red' color with a dashed linetype.The thickness of the lines were mapped based on the weight value Fig. 7 qRT-PCR verification of DEGs.Seven genes were randomly chosen for qRT-PCR validation.Data was shown as mean ± SD (standard deviation).The bars represented the standard deviation (n = 3).The actin was used as an internal control to normalize the expression data.Note: In A-G the left of the double axis represented RNA-seq, the right of the double axis represented the qRT-PCR; In H, correlation between the relative expression of transcriptome and qRT-PCR Plant varieties with drought resistance have armed with higher antioxidant enzyme activities [41].Drought induced an increase in TAOC which is proportional to the duration and intensity of water deprivation [27].The results showed that the TAOC of P. pygmaea was significantly higher than that of the other two pines (Fig. 1A), which may have stronger drought tolerance.After rehydration, the TAOC of the three species remained at a relatively high level.This maintenance of high TAOC after rehydration is similar to other plants [42].The SP of P. elliottii exhibited an initial decrease followed by an increase, with levels significantly higher than those of P. massoniana and P. pygmaea at several points between 14 and 35 days of drought stress (Fig. 1B).Above-ground tissues of drought trees showed increases in both sucrose and starch [43].NSC are major substrates for plant metabolism and have been implicated in mediating drought-induced tree mortality [44].NSC are important "buffers" to maintain plant functions under drought conditions, the seedlings under the severe and moderate drought conditions tended to invest newly fixed C to starch storage and hydraulic repair instead of growth [45].Once drought ceased, C allocation to storage was still prioritized at the expense of growth [46].Severe drought treatment largely reduced photosynthetic carbon assimilation [47].Plants may have allocated more NSC from source organs to roots [48].After rehydration, the content of starch and NSC in P. pygmaea was higher than that of the other species, indicating a strong recovery ability (Fig. 1B).The occurrence of drought stress can affect the decrease of photosynthetic parameters and chlorophyll fluorescence [27].The SS content, starch content, and NSC (Fig. 1C-E) of all three pines displayed a "W" pattern.Drought stress lasted for 15 days caused obvious photosynthesis inhibition in Atractylodes lancea [27].In general, drought significantly reduce Chla and Chlb content [49].During drought, the content of Chla and Chlb decreases continuously (Fig. 1H&I).The activities of CAT and POD increased under drought stress [50], the activities of POD and SOD increased significantly compared with the control [51].After rehydration, P. pygmaea's POD activity was higher than the other two species (Fig. 1F), and MDA content was lower than P. elliottii and P. massoniana at Re7D (Fig. 1G).Most antioxidant enzyme indicators still fail to reach the control level after 7 days of rehydration, which may indicate that the 7 day period is insufficient to offset the impact of drought.Other study have shown that, the recovery of P. massoniana MAD needs to continue to decrease to the control level after 15 days of rehydration [52].Long term drought still requires a period of recovery, this recovery time may be determined by the characteristics of species.
GSEA showed that "Linoleic acid metabolism", "Alpha-Linolenic acid metabolism" were upregulated in the Pp35D (compared to the Pp0D)."Cytoskeleton proteins", "Glycosaminoglycan binding proteins", "Photosynthesis-antenna proteins" were downgulated in the Pp35D (Fig. 4I).Alpha-Linolenic acid significantly increases under drought stress [61].Compared with fully watered P. taeda seedlings, levels of some unsaturated fatty acids significantly reduced in needles of under water stress, such as palmitelaidic acid, erucic acid, and alpha-linolenic acid [62].The needles of P. massoniana seedlings may respond to drought mainly through regulating ABA and JA hormone-related pathways [63].ABA can accumulate 10-to 30-fold in plants under drought stress relative to unstressed conditions [64].ABA signaling pathway is an important signal, induces stomatal closure and inhibition of transpiration under drought stress [65,66].ABA-responsive genes were obviously upregulated at drought treatment for 7 days [67].The GO term of "abscisic acid-activated signaling pathway" was enriched in two drought term (Pp28D, Pp35D) compared with PpCK (Fig. 4J).Two SNRK2 (sucrose non-fermenting 1-related protein kinase 2) genes, five PYL genes, eight PP2C genes, seven ACAA1 genes and three ABF genes were enriched in the core enrichment of the between Pp35D and PpCK.PYR/PYLs are ABA receptors functioning at the apex of a negative regulatory pathway that controls ABA signaling by inhibiting PP2Cs [68].PYLs/ PP2Cs heterodimer blocks substrates binding to PP2Cs, and thus releases SnRK2s whose kinase activity is formerly inhibited by PP2Cs through physical interaction and dephosphorylation, finally allows the plants to adapt to the stress conditions [69].TaABFs-regulated PYL gene (TaPYL1-1B) exhibited higher ABA sensitivity, photosynthetic capacity and water-use efficiency, all contributed to higher drought tolerance than that of wild-type [70].Alpha-linolenic acid mostly participated alleviated drought stress through JA signaling [71].As a precursor to the synthesis of JA, α-linolenic acid in lipid metabolism is related to drought stress [71].Different strains of tobacco seedlings were enriched in functions about alpha-linolenic acid, and arginine and proline metabolisms under drought [72].In the alpha-linolenic acid pathway, two alcohol dehydrogenase (ADH1) genes, five Lipoxygenase (LOX) genes were enriched.LOX play essential roles in responding to biotic and abiotic stresses through oxidizing various lipids [73].Overexpression of oriental melon CmLOX13 improved drought tolerance in Arabidopsis [74].Compared with the Pp35D, "Diterpenoid biosynthesis", "Prodigiosin biosynthesis", "Nicotinate and nicotinamide metabolism", "Terpenoid backbone biosynthesis", "Photosynthesis", "Fatty acid biosynthesis", "Fatty acid biosynthesis", "Lipid biosynthesis proteins", "Glycosyltransferases" were upregulated in the PpRe7D (Fig. 4K).These findings suggest that P. pygmaea may respond to drought by enhancing metabolic processes such as ABA signaling pathway, alpha-linolenic acid, which enhances overall antioxidant capacity.
The results of WGCNA showed that bisque4 module was positively associated with POD (Fig. 5A).In the bisque4 module, "Flavonoid biosynthesis" (ko00941) and other pathways were enriched.Three GST gene were enriched in the bisque4 module.Overexpression of MsGSTU8 of Medicago sativa in Nicotiana tabacum resulted in maintained chlorophyll content, increased antioxidant enzyme activity, and soluble sugar levels after salt-alkali treatments [75].These findings suggest that GST genes play an important role in scavenging reactive oxygen species and enhancing antioxidant enzyme activity.The brown module was negative correlation with TAOC.The genes in the brown module were significantly inhibited in expression after drought."Cytoskeleton proteins" (ko04812), "Photosynthesis-antenna proteins" (ko00196), "Phenylpropanoid biosynthesis" (ko00940), "Porphyrin and chlorophyll metabolism" (ko00860) and others were enriched.This indicates that drought can inhibit the synthesis of genes involved in photosynthesis, chlorophyll, and phenylpropanoid biosynthesis.The plum2 was positively associated with TAOC.Four CAT genes, protein transport SEC23, LEC14B homolog were enriched in the plum2 module.The cultivars Coral Orange, Tendersweet and Solar Yellow were tolerant to drought stress, which was supported by their higher transcript levels of CAT , SOD genes [76].The drought tolerant cultivars have higher accumulated of CAT gene [77].LEC14B (JK340642) had significantly higher level of expression in drought-tolerant genotype (Tifway) than in drought-sensitive genotype (C299) at 10 d of drought, but no significant differences were detected between two genotypes at 5 of drought [78].These suggests that CAT, LEC14B genes may enhance antioxidant capacity and suppress ROS.Sec23/Sec24 are important cytosolic components of vesicle coat protein II (COPII) [79].SEC23 is responsible for transporting newly synthesized proteins and lipids from the endoplasmic reticulum (ER) to the Golgi apparatus in cells, Steptoe SEC23 (A0A287NVF2) was down regulated under N deficiency, drought stress or both stress [80].In our study, two SEC23 (D1_tran-script_5780, D1_transcript_28060), two SEC13, two SEC24, three SEC31 were enriched in the brown modules, showing a relatively low expression trend in the later stages of drought (Fig. 5D).But not all SEC genes expression were inhibited, the SEC23 (D1_transcript_53479) was one of the plum2 module and one of the Profile 39 (0, 1, 2, 3, 4, 5, 6) genes, showed a upregulated trends under drought stress, meaning SEC23 in P. pygmaea play important roles in regulating both trafficking and plasma membrane stability.
Considering the significance of the selected module hub genes, we will proceed to validate the functions of the related genes through gene cloning and overexpression.Furthermore, it is worth noting whether the drought-tolerant type of P. pygmaea, when used as a rootstock for grafting with P. massoniana scions, leading to changes in phenotypes such as scion height and resistance in the grafted plants.This study will provide a valuable direction for seedling nursery establishment, rootstock breeding, and reduction of artificial production costs in conifer orchards.

Conclusions
Different species exhibit inconsistent drought tolerance.Physiological assessments of drought and rehydration responses in three pine species (P.pygmaea, P. elliottii and P. massoniana) indicate that P. pygmaea exhibited a stronger drought response and a strong recovery ability under prolonged drought conditions.High-quality transcriptome of P. pygmaea was obtained using PacBio SMRT sequencing.Comparative transcriptome analysis indicates that drought suppresses pathways related to photosynthesis and chlorophyll synthesis metabolism.P. pygmaea may respond to drought by enhancing metabolic processes such as ABA signaling pathway, alpha-linolenic acid, which enhances overall antioxidant capacity.WGCNA revealed GST, CAT , LEC14B, SEC23 were associated with antioxidant enzyme activity and TAOC.qRT-PCR analysis correlated well with the data obtained by RNAseq.This study provides valuable insights into the drought response and regulatory processes of P. pygmaea.

Plant material and drought treatment
P. pygmaea exhibit a clear absence of apical dominance with an inconspicuous main trunk.The base is multibranched, clumped, with a height range of 40-50 cm to 1-2 m.Leaves grow short and upright while its seed cones cluster together [3].Tracheid length of P. pygmaea is shortened, resin canals increase and more or less concentrate in the late wood [7].P. elliottii trees grow to 30 m tall.Through early seed resource investigation and collection, the seeds of half-sib progeny of P. pygmaea, P. massoniana and P. elliottii were obtained from the highquality half-sibling progeny of Ma'anshan Seed Garden (26°16' N and 107°31' E), Duyun City, Guizhou Province, China.The seeds were planted in the greenhouse, using a humus: yellow loam soil mixture (1:3) as the soil medium.In order to obtain a detailed genetic background of P. pygmaea, different parts of a 5-year P. pygmaea were collected for PacBio SMRT sequencing, including the roots, needles, ovulate strobilus, and staminate strobilus.We selected three types of annual pine trees as the experimental subjects.Prior to the drought treatment, the three types of pine trees were irrigated normally.Sampling was conducted between 8:00-9:00 in the morning, mature needles were collected from the stem tips every seven days from new plants, labeled as 0D, 7D, 14D, 21D, and 28D during the drought treatment.A total of 35 days between the cessation of irrigation and the end of the severe drought.When the drought has lasted for 35 days, some seedlings began to die, and the drought treatment is stopped (labeled as 35D).After rehydration, the seedlings were left to grow for seven days.At the end of this period, samples were collected and labeled as Re7d.Three biological replicates were used for each treatment.

Determination of physiological indexes
The TAOC of needles was determined using 2,2'-azinobis[3-ethylbenzothiazoline-6-sulfonic acid] diammonium salt (ABTS •+ ) method [81].The frozen samples (about 2.0 g) were homogenized with 15 mL of 50% methanol and centrifuged at 10,000 rpm for 20 min at 4 °C, the supernatant was collected for ABTS analysis.SP content was measured by Coomassie Brilliant Blue G-250 staining [82].Fresh needles (0.1 g) were cut and fully ground in phosphate buffer.Then added 5 mL of Coomassie Brilliant Blue G-250 solution.The mixture was then fully mixed and incubated for 2 min at room temperature.The absorbance at 595 nm was measured.SS content was determined using phenol method [83].Fresh needles (0.5 g) was extracted in 10 mL of ddH 2 O for 20 min at 90 °C and centrifuged at 4000 rpm for 5 min.A combination of 1 mL of sample extract, 0.5 mL of anthroneethylacetate reagent, and 5 mL of concentrated sulfuric acid was homogenized and boiled for 5 min, and then rapidly cooled.The absorbance of the mixtures was detected at 630 nm.Starch content was measured by the anthrone colorimetric [84].NSC was defined as the sum of soluble sugar and starch concentrations [85].POD enzyme activity was determined using guaiacol method.The reaction mixture was as follows, 50 mmol/L Tris-HCl buffer (pH 7.0) containing 0.1 mmol/L EDTA, 10 mmol/L guaiacol, 5 mmol/L H 2 O 2 and 100 μL enzyme extract.After the enzyme solution was added to the reaction system, it was immediately incubated at 34 °C in a water bath for 3 min.The changes in absorbance values were measured at 470 nm.One unit of POD activity was expressed as units (μmol guaiacol decomposed per minute) per mg of fresh weight (FW).MDA content was determined using thiobarbituric acid (TBA) method [86].Each fresh sample (about 0.5 g) was homogenized with 4 ml of 0.1% 0.1% (w/v) trichloroacetic acid (TCA) in ice bath and centrifuged at 4000 rpm for 10 min at 4 °C, for every 1 ml aliquot, 1 ml of 20% (w/v) TCA comprised of 0.5% (w/v) TBA was added.The mixture was heated at 95 °C for 30 min and then cooled immediately before centrifugation again.The absorbance of the supernatants was detected at 532, 450, and 600 nm.Fresh needles were cut into pieces and merged in 80% (v/v) acetone in the dark at room for 24 h to extract chlorophyll.Chla, Chlb were calculated by ultraviolet-visible spectrophotometer at 663, 646 and 470 nm [87].Three biological replicates were used for the above experiments.

PacBio library construction,sequencing, and annotation
Total RNA was extracted according to the instruction manual for Trizol reagent (Invitrogen, Carlsbad, CA, USA).RNA integrity was assessed by agarose gel electrophoresis, RNA purity (OD260/280 and OD260/230), concentration and RNA integrity number (RIN), 28S/18S were detected using a NanoDrop 2000 spectrophotometer and an Agilent 2100 (Agilent Technologies, Santa Clara, California, USA).The mRNA was enriched with Oligo (dT) magnetic beads.After all RNA samples were qualified, the total RNA from different tissues was mixed at equal ratios.After the library was qualified, it was sequenced base on the PacBio platform (Menlo Park, CA, USA).The raw bam file had been deposited in the NCBI SRA database (BioProject acc.PRJNA994856).The raw seqencing subreads were filtered using SMRTLink v8.0 (https:// www.pacb.com/ suppo rt/ softw are-downl oads/) with default parameters.The CCS was obtained by merging the subreads from the same polymerase reads, CCS were obtained by self-correction.FLNC sequences refer to a kind of CCS that contains both 5' and 3' primers and poly-A tail, without chimeric reads.After clustering the FLNC reads using the IsoSeq module, the consistent sequences in each cluster were further corrected.Finally, high-quality full-length transcripts (HQ) with an accuracy greater than 99% and low-quality (LQ) transcripts were obtained.Using CD-HIT software [88] to remove redundant sequences in transcripts and cluster the transcript sequences.BUSCO v5.3.2 [89] was used to assess the assembly quality, the BUSCO dataset of 1603 predefined single-copy genes dedicated to gymnosperms [90] were download and used.The transcripts were annotated using the NCBI non-redundant protein (NR) database, the Swiss-Prot protein (http:// www.expasy.org/ sprot/), the Eukaryotic Orthologous Groups (KOG) database (http:// www.ncbi.nlm.nih.gov/ COG/), the Pfam database (http:// pfam.sanger.ac.uk/), Kyoto Encyclopedia of Genes and Genomes (KEGG) (http:// www.genome.jp/ kegg, cutoff e-value ≤ 1e-5) pathways, and the eggNOG database (http:// eggno g5.embl.de/#/ app/ home) [91].The SSRs were identified using the MISA software [92].To identify long non-coding RNA (LncRNA), the transcripts with coding potential were filtered using four computational approaches, including the Coding Potential Assessment Tool (CPAT), Coding-Non-Coding Index (CNCI), Coding Potential Calculator (CPC), and Pfam protein structure domain analysis (Pfam).

RNA-seq library preparation, sequencing
The needles samples of different stages of P. pygmaea's response to drought and rehydration were collected.The total RNA was isolated using the Trizol Reagent (Invitrogen, Carlsbad, CA, USA).RNA purity, concentration and integrity were assessed using agarose gel electrophoresis, NanoDrop 1000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA) and Agilent 2100 (Agilent Technologies, Santa Clara, California, USA).The mRNA was enriched with Oligo (dT) magnetic beads.The mRNA was added to fragmentation buffer and cut into short fragments.Using mRNA as a template, cDNA was reverse-transcribed using six-base random primers.The double-stranded cDNA samples were purified, endrepaired, added with poly(A) tails, and then ligated to sequencing adapters to create cDNA libraries, qualified libraries were sequenced using the NovaSeq 6000 system with paired-end reads.The raw reads generated from Illumina sequencing were deposited in the NCBI SRA database (BioProject acc.PRJNA972125).

Expression quantification, identification of DEGs
Fastp v0.20.1 was used for for quality control [93].The non-redundant transcripts which obtained from Pacbio SMRT sequencing was used as background transcripts.Bowtie2 v2.4.3 [94] was used to align the sequenced paired-end transcriptomic data.The gene expression level was determined according to the fragments per kilobase of transcript per million mapped reads (FPKM) by the RSEM v1.2.12 [95].The DEGs were identified based on the counts using DESeq2 v1.40.1 R package [96].DEGs screening threshold was set to adjusted p-value cutoff < 0.05 and |foldchange|> 2. Gene Ontology (GO) and KEGG enrichment analysis were performed by the clusterProfiler v4.8.1 R package [97].STEM software [98] was used to perform the trend clustering analysis.GSEA based on GO terms and KEGG pathways were conducted by the clusterProfiler, gene sets were considered statistically significant at an p.adjust of 0.05.Gene co-expression networks were constructed using the WGCNA v1.72.1 R package [99].All the DEGs were filtered with at least one sample FPKM value greater than 5.

Quantitative real-time PCR analysis
In order to verify the accuracy of transcriptome, seven genes were randomly selected for qRT-PCR.The Primer Premier 5.0 software was used to design the primers (Table S1).The first strand of the cDNA was synthesized from total RNA.qRT-PCR was performed on a realtime CFX96 Touch Real-Time PCR System (Bio-Rad, USA).Each 50 μL reaction included: 25 μL of 2X SYBR ® Green PCR Master Mix, 1μL of Forward /Reverse primer (10 μM), 2.0 μL of cDNA template, and 21 μL of RNasefree water.The PCR reaction conditions were as follows: preheating at 95 °C for 30 s, 40 cycles of heat denaturation at 95 °C for 5 s, annealing at 60 °C for 34 s.The actin gene was used as a reference gene.The gene relative expression levels were calculated according to the 2 −ΔΔCt method [100].

Statistical analysis
The variance of the data was analyzed using the R v4.2.3 software, and the significance threshold was set at p < 0.05.LSD was used at a probability level of 0.05 to verify the significance between species.PCA was performed using the FactoMineR v2.8 R package [101].

Fig. 1
Fig. 1 The physiological changes of three pines under drought and rehydration.A Total antioxidant capacity (TAOC); B Soluble proteins (SP); C Soluble sugars (SS); D Starch; E Non-structural carbohydrate (NSC); F Peroxidase (POD) activity; G Malondialdehyde (MDA) content; H Chlorophyll a (Chla); I Chlorophyll b (Chlb); J The relative change rate of indicators compared to CK at different time points (ratio = (Value (point) -Value (CK) )/Value (CK) ); K PCA analysis.Note: In (A-I), 7-35D represented drought treatment after 7 D, 14D, 21D, 28D, and 35D.Re7D represented 7 D after rehydration; error bars indicated the standard deviation, the number of seedlings measured for each species was n = 3; the left Y-axis represented value of physiological traitors, "CK" group of each species was set as a reference to compare the mean by the Student's t-Test, "ns" means no difference, "*" P ≤ 0.05;"**" P ≤ 0.01;"***" P ≤ 0.001;"****" P ≤ 0.0001.A least significant difference test (LSD) was used at a probability level of 0.05 to verify the significance between species at sample point

Fig. 3
Fig. 3 Transcriptome data analysis and differential expression gene quantity analysis.A PCA; B Summary of the identified DEGs; C Venn diagram of DEGs between PpCK and other groups; D Venn diagram of DEGs between Pp35D and other groups; E The clustered and heatmap analysis of DEGs